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I A NEW  APPROACH  TO  THE  ANALYSIS  OF 

STOCHASTIC  LANCHESTER  PROCESSES:  I.  TIME  EVOLUTION 

Abstract 

A new  approach  to  the  study  of  stochastic  Lanchester  pro- 
cesses based  on  diffusion  approximations  is  presented.  The 
distribution  of  the  two  force  levels  over  time  is  shown  to  be 
well  approximated  by  a nonstationary  bivariate  Gaussian  diffu- 
sion process  with  specified  mean  and  covariance  structure.  The 
approximation  is  based  on  an  asymptotic  analysis  which  assumes 
the  initial  force  levels  are  large.  Numerical  studies  are 
presented,  however,  which  show  surprising  accuracy  for  force 
levels  as  small  as  30. 

A wide  variety  of  attrition  stiructures  are  discussed  in- 
cluding the  linear  and  square  law  cases,  Helmbold' s general 
attrition  structure,  Karr' s engagement  model,  and  heterogeneous 
models.  The  development  of  tractable  mathematical  expressions 
for  the  time  evolution  of  complicated  Lanchester- type  attrition 
processes  makes  possible  the  Introduction  and  auialysis  of  de- 
cision theoretic  aspects  to  the  problem  such  as  force  level 
decisions,  combat  tactics,  reinforcement  decisions,  and  the 
value  of  information  about  the  opponent' s strengths,  weaknesses, 
and  strategies.^ 


1.  Introduction 

The  mathematical  theory  of  the  attrition  of  forces  engaged 
in  combat  was  initiated  by  Lanchester  [12]  in  1914.  In  his 
pioneering  work,  the  attrition  process  was  formulated  as  being 
deterministic  and  was  described  by  a set  of  ordinary  differential 
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equations.  Many  advsLnces  have  been  made  beyond  this  work.  These 
advances  include  such  complicaticns  as  time  dependencies,  complex 
attrition  structures,  force  inhomogeneities,  and  aspects  of  control. 
The  reader  should  consult  the  tutorial  by  Taylor  [14]  for  a survey 
of  this  work  and  a complete  bibliography.  It  is  fair  to  say 
that  most  of  the  advances  have  come  within  the  deterministic 
point  of  view,  while  the  stochastic  theory  has  been  much  slower 
to  develop. 

Within  a stochastic  framework,  some  of  the  more  importaint 
contributions  are  due  to  Clark  [4],  Karr  [10],  Weiss  [15],  Brooks 
[3],  and  Helmbold  [9].  Both  [15]  and  [5]  present  the  idea  that 
under  certain  conditions  a stochastic  analysis  is  not  needed, 
because  the  results  differ  little  from  that  of  a deterministic 
analysis.  Karr  [10]  presents  a methodology  for  developing 
Markov  chain  models  from  explicit  combat  assumptions.  Using 
this  methodology,  a large  number  of  models  are  presented;  how- 
ever, none  is  analyzed.  An  immediate  problem  encountered  with 
a stochastic  analysis  is  that  the  interesting  models  which 
arise  usually  cannot  be  solved,  at  least  exactly.  In  fact,  only 
Clark  [4]  presents  closed  form  expressions  for  the  "linear"  law, 
and  these  are  very  complex  euid  involve  a large  number  of  expo- 
nential terms. 

In  spite  of  the  lack  of  development  of  such  a theory  and 
the  mathematical  intractabilities  encountered,  there  are  two 
very  important  reasons  for  adopting  a stochastic  approach. 

First,  we  take  issue  with  the  notion  of  "stochastic  determinism" 
presented  in  [5]  and  [15]*  The  claim  has  been  made  that  for 
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reasonably  large  Initial  force  levels,  quantities  such  as  the 
probability  of  victory  become  nearly  0 or  1 and  reflect  the  re- 
sults of  the  deterministic  analysis.  Nevertheless,  other 
quantities  are  important  and  exhibit  large  random  fluxuations. 

For  example,  the  force  levels  themselves  plotted  over  time 
will  exhibit  enormous  variability.  This  variability  must  be 
characterized,  especially  if  one  wishes  to  introduce  elements 
of  tactics  or  strategies  which  might  depend  on  tire  actual  evo- 
lution of  the  combat.  This  variability  can  be  seen  in  the 
results  presented  in  section  4. 

The  second  reason  for  using  a stochastic  analysis  is  that 
it  allows  one  to  introduce  decision  theoretic  concep  s to  the 
study  of  Lanchester  processes  and  hence  to  introduce  a new  role 
for  this  area  of  operations  research.  If  Lanchester  theory  is 
to  have  practical  value,  it  should  allow  one  to  address  problems 
such  as  the  value  of  information  about  the  opponent's  force 
levels  and  strategies,  what  type  of  battle  to  fight,  when  to 
call  in  reinforcements,  and  when  to  avoid  conflict  altogether. 

Some  of  these  questions  are  addressed  in  Schreiber  [13] • The 
answers  to  these  questions  require  the  introduction  of  a decision 
theoretic  structure  to  the  problem,  and  they  require  a stochastic 
approach.  Of  course,  a stochastic  approach  is  not  enough.  One 
must  be  able  to  obtain  tractable  mathematical  expressions  which 
can  then  be  used  in  a decision  analysis. 

In  this  paper,  a new  approach  to  stochastic  Lanchester  theory 
is  Introduced.  It  is  approximate  in  nature  and  is  based  on  the 
method  of  diffusion  approximations.  The  reader  should  consult 
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Gaver  and  Lehoczky  [5]  and  [6]  for  other  applications  of  this 
methodology.  The  method  is  asymptotic  and  provides  results 
ass-uning  the  initial  force  levels  are  large.  Nevertheless, 
numerical  results  will  be  presented  which  indicate  that  the 
method  offers  excellent  accuracy  for  initial  force  levels  as 
small  as  JO.  This  accuracy  when  coupled  with  the  analytic 
tractability  of  the  expressions  derived  make  the  method  ideal 
for  studying  random  fluxuations  as  a function  of  the  relevant 
parameters  and  tactical  aspects  of  Lanchester  theory. 

The  diffusion  approximation  approach  is  Introduced  in  sec- 
tion 2 and  is  applied  to  the  "square"  law  model.  Section  J 
contains  an  analysis  of  a number  of  different  attrition  struc- 
tures including  nonlinear  structures  and  heterogeneous  models. 
Section  4 presents  numerical  results  which  assess  the  accuracy 
of  the  approximation.  The  application  of  this  methodology  to 
solving  tactical  decision  problems  is  the  subject  of  another 
report. 
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2.  The  Diffusion  Approximation  Approach 

In  this  section,  Markov  chain  models  of  Lanchester  processes 
are  introduced  and  analyzed  using  diffusion  approximations. 

Only  the  simple  "square  law"  studied  in  this  section  in  order 
to  illustrate  the  methodology.  Complex  models  including  non- 
linear attrition  structures  and  heterogeneous  models  will  be 
examined  in  a subsequent  section.  Throughout  this  section,  we 
assume  there  are  two  opponents  each  with  homogeneous  forces 
which  undergo  attrition.  If  x(t)  and  y(t)  denote  the  two 
force  levels  as  a function  of  time,  then  classical  Lanchester 
theory  would  model  the  functions  x and  y as  deterministic 
and  satisfying  a system  of  differential  equations 

x'(t)  = -f(x(t),  y(t))  , x(0)  = x^,  y(0)  = y- 

° ° (2.1) 

y (t)  = -g(x(t),  y(t)). 

The  functions  f and  g are  chosen  to  model  the  type  of 

conflict  involved.  For  example,  if  f(x,  y)  = axy  and  g(x,  y)  = 

bxy,  then  the  model  is  called  the  "linear  law"  (since  ^(xq- 

x(t))  = a(yQ-y(t)))  and  was  Introduced  by  Lanchester  as  a model 

of  "ajccient"  combat.  If  f(x,  y)  = ay  and  g(x,  y)  = bx, 

2 2 

then  the  model  is  called  the  "square  law"(since  b(xQ-x  (t))  = 

^(Vq  - y (t)))  and  was  introduced  as  a model  of  "moaern"  combat. 
Many  functions  f and  g have  been  proposed  to  model  various 
aspects  of  combat.  The  reader  should  consult  the  tutorial  by 
Taylor  [l4]  for  a survey  of  such  models  and  a large  bibliography. 

The  functions  f and  g can  also  depend  explicitly  on  time. 
We  do  not  allow  for  this  possibility  in  this  paper,  but  it  should 
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be  pointed  out  that  the  method  of  diffusion  approximations  can 
easily  accomodate  it. 

Stochastic  models  of  attrition  processes  are  typically 
formulated  as  continuous  time  Markov  chains  with  state  space 
S = [ ( i,  j)|i>0,  J>0}.  The  deterministic  system  given  by 
(2.1)  has  a natural  stochastic  analogue.  Specifically,  let 
X(t)  and  Y(t)  represent  the  two  force  levels  at  time  t. 

We  assxime  C(X(t),  Y(t)),  t>0]  is  a Markov  chain  with  transi- 
tions and  transition  rates  given  by 


t t + dt 

Rate 

(i,  j)  (i-i,  j) 

f(i,  J)dt  + o(dt) 

(i,  J -1) 

g(i,  J)dt  + o(dt) 

(2.2) 

(1,  J)  l-(f(l,  J)+g(i,  J))dt  + o(dt) 

all  others  o(dt). 

If  we  define  dX(t)  = X(t  + dt)  - X(t)  and  dY(t)  = Y(t  + dt)  - 
Y(t),  then  (dX(t),  dY(t))  is  a random  vector  with  mean  (-f(X(t), 
Y(t))dt,  -g(X(t),  Y(t))dt)  and  covariance  matrix 

f(X(t),  Y(t))dt  0 \ 

0 g(X(t),  Y(t))dt  / 

to  terms  of  order  dt.  For  the  deterministic  process  specified 
by  (2.1),  the  vector  (dx(t),  dy(t))  is  (-f(x(t),  y(t)), 

-g(x(t),  y(t))),  so  the  two  processes  are  adapted  through  the 
mean.  The  Markov  chain  model  is  a generalization  of  the  deter- 
ministic model  in  that  (2.1)  can  be  recovered  in  the  mean  but 
it  also  takes  Into  account  the  fluxuatlons  encountered  in  real 
attrition  situations. 
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The  Markov  chain  model  specified  by  (2.2)  brings  with  it 
immediate  complications.  To  specify  the  distribution  of  (X(t), 
Y(t))  exactly,  one  needs  to  solve  the  Kolmogorov  forward  equations: 

= f(m+l,  n)P(t,  m+1,  n)  + g(m,  n+l)P(t,  m,  n+1) 
dt  (2.3) 

- (f(m,  n)  + g(m,  n))P(t,  m,  n), 

where  0<m<XQ,  0 < < Yq 

P(t,  m,  n)  = P((X(t),  Y(t))  = (m,n)l(x(0),  Y(0))  = (xq,  y^)). 

Equations  (2.3)  are  essentially  impossible  to  solve  in  any 
convenient  closed  form.  Only  Clark  [4]  has  developed  exact 
expressions  for  P(t,  m,  n);  however,  they  were  only  for  the 
linear  law  case,  and  the  results  involve  a double  sum  of  expo- 
nentials. Such  exact  expressions  are  of  little  value  in  pro- 
viding information  about  the  behavior  of  the  process  as  a function 
of  the  initial  values  and  the  attrition  coefficients.  If  one 
wishes  to  adopt  a Markov  chain  approach,  an  alternate  method  of 
analysis  must  be  developed. 

The  method  of  diffusion  approximations  provides  such  an 
approach  when  the  initial  force  levels  Xq  and  y^  are  large. 

The  state  space  S is  replaced  by  a continuous  state  space,  and 
the  sample  paths  become  continuous  functions  rather  than  step 
functions.  If  X(t)  and  Y(t)  are  large,  then  over  any  short 
time  period  many  transitions  will  take  place.  Since  the  total 
change  in  the  X and  Y processes  will  be  the  sum  of  many  unit 
steps,  the  total  change  will  be,  by  virtue  of  the  central  limit 
theorem,  approximately  normally  distributed.  The  first  step  in 
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the  approximation  is  to  replace  (2.1)  and  (2.2)  by  a diffusion 
process  defined  by  the  Ito  stochastic  differential  equations 
(s.d.e.  in  the  sequel) 

dX(t)  = -f(X(t),  Y(t))dt  + (f(X(t),  Y(t)))^/^dW.  (t) 

(2.4) 

dY(t)  = -g(X(t),  Y(t))dt  + (g(X(t),  Y(t)))^/^dW2(t) 

with  X(0)  = Xq,  Y(0)  = Yq  a.s.,  and  where  W^(t)  and  W^(t) 
are  independent  standard  Wiener  processes. 

The  increments  of  {(X(t),  Y(tl),  t>0]  are  normally 
distributed.  The  s.d.e.  (2.4)  corresponds  to  a partial  differ- 
ential equation,  the  Kolmogorov  forward  equation,  which  is  similar 
to  (2.3)  and  no  easier  to  solve.  Nevertheless,  an  approximate 
analysis  can  be  carried  out  if  Xq  and  y^  are  large.  Speci- 
fically, we  let  X(0)  = Nu(0)  and  Y(0)  = Nv(0).  We  introduce, 
in  the  spirit  of  the  central  limit  theorem,  U^j(t)  = (X(t)  - 
Nu(t))/N^/^  and  Vjj(t)  = (Y^(t)  - Nv(t))/N^/^.  The  functions 
u(t)  and  v(t)  are  deterministic.  The  process  [(U^(t),  Vjj(t)), 
t>0}  will  also  be  a diffusion  process,  and  the  system  of  s.d.e. 's 
governing  it  can  be  found  using  Ito's  Lemma  (Arnold  [1],  p.90  or 
Gihkman  and  Skorokhod  [7]>  P*27)  to  be 

dUjj(t)  = -N^/^(u'(t)  + f(Nu(t)  + N^/^Ujj(t),  Nv(t)  + N^/^Vjj(t)))dt 
+ (f(Nu(t)  + N^/^u^(t),  Nv(t)  + N^/^Vjj(t)))^/2dW^(t) 

(2.5) 

dVjj(t)  = -N^/2(v'(t)  + g(Nu(t)  + N^/^Ujj(t),  Nv(t)  + N^/^Vj^(t)))dt 
+ (g(Nu(t)  + N^/^Uj^(t),  Nv(t)  + N^'^\(t)))^^^dW^(t) 


with  Uj,(0)  » Vjj(O)  =0  a.  s. 
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We  wish  to  study  (2.5)  assuming  N->oo;  however,  to  proceed 
further  we  must  make  specific  assumptions  about  f and  g.  In 
this  section,  we  assiame  f(x,  y)  = ay  snd  g(x,  y)  = bx,  the 
quadratic  law  case.  In  general  we  assume  f(X(t),  Y(t))  and 
g(X(t),  Y(t))  are  of  order  N.  This  allows  for  an  asymptotic 
analysis  of  (2.5)  based  on  Theorem  K of  Barbour  [2]  or  on 
Kurtz  [11]. 

In  the  specific  quadratic  case,  equations  (2.5)  become 

dUj^(t)  = -N^/^(u'(t)  + av(t))dt  - aYjj(t)dt  + 0(N~^/^)dt 

+ (av(t)  + 0(N"^/^) )^/^dW.  (t) 

^ (2.6) 

d\(t)  = -N^/^(v'(t)  + bu(t))dt  - bX^(t)dt  + 0(N”^^^)dt 
+ (bu(t)  + 0(N‘^/^))^/^dW2(t) 


with  Uj^(O)  = V^(0)  = 0 a.s. 

Clearly,  the  coefficients  of  the  terms  must  be 

identically  0 for  {(Ujj(t),  Vjj(t)),  1 2’0}  to  converge  to  a 
sensible  limit  as  t-^oo.  Indeed,  Theorem  K of  [2]  shows  that 
{(Ujj(t),  Vj^(t)),  t>0}  converges  weakly  to  a limiting  diffusion 
process  {(U(t),  V(t)),  t^O]  satisfying  the  system  of  Ito 
s.d.e.'s  given  in  matrix  form  by 


/dU(t)\ 

\dV(t)] 


U(t)  dt  + |(av(t))^/2 

\v(t)  I \ 0 


or 


dU(t)\  ^ 

At 

dt  + B. 

f dW^(t)  ’ 

\dvit)  1 " 

t 

dW2(t)j 

0 \ [dW^(t)\ 

(bu(t))^/2  I ^dW2(t)| 

(2.7) 


with  U(0) 


V(0)  = 0 a.s.. 
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provided  u and  v satisfy  the  system  of  ordinary  differential 
equations 

u'(t)  = -av(t)  u(0)  = Xq,  v(0)  = y^ 

V' ( t)  = -bu( t)  (2.8) 

Of  course  (2.8)  are  Just  the  ordinary  deterministic  Lan- 
chester  equations.  The  diffusion  approximation  method  consists 
of  replacing  (X(t),  Y(t))  which  is  identically  equal  to 
N(u(t),  v(t))  + N^/^(Ujj(t),  Vjj(t))  by  N(u(t),  v(t))  + N^/^(U(t), 
V(t)),  where  the  latter  is  given  by  (2.7)  and  (2.8)  and  is  far 
easier  to  characterize  then  the  former.  This  approximation  is 
based  on  exactly  the  sajne  idea  as  the  central  limit  theorem 
which  would  approximate  a binomial  (N,  p)  distribution  by 
Np  + Normal  (O,  p(l-p)).  The  first  two  moments  have  been 

fit  exactly,  and  probability  statements  are  exactly  correct  in 
the  limit  and  very  accurate  for  even  samll  values  of  N pro- 
vided p is  not  too  extreme.  Numerical  comparisons  will  be 
given  later  which  show  that  this  method  provides  excellent  a 
accuracy  for  values  of  N as  small  as  ^0. 

The  method  is  useful  because  the  (U(t),  V(t))  process 
can  be  explicitly  characterized.  The  system  of  s.d.e.'s  (2.7) 
is  linear,  and  much  is  known  about  such  systems  (for  example  [1], 
chapter  8).  Specifically,  E(U(t),  V(t))  = (0,0)  for  all  t 
and  (U(t),  V(t))  will  have  a bivariate  normal  distribution 
for  all  t.  The  process  is  thus  a nonstationary,  bivariate, 
Gaussian  diffusion  with  mean  0 and  covariance  matrix  = 
E((U(t),  V(  t)  )'^(U(  t) , V(t)))  which  is  the  unique  non-negative 
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definite  solution  of  the  matrix  Ricatti  equation 


It  ~ -t^t  ~t~t  -t2t^  ^0  - 2 


Equations  (2.9)  can  also  be  written 


(2.9) 


Oll(t) 

o^2(t) 

^) 


0 -2a  0 


0 -a 


0 -2b  0 


'^12^^) 

c^22^ 


av(  t) 


bu(  t) 


or  o(t)  = c^o(t)  + 

with  Oj_j(0)  = 0,  a^^(t)  = Var(U(t)),  ^ Var(V(t)),  and 

a^2(t)  = Cov(U(t),  V(t)). 

Equations  (2.8)  and  (2.10)  can  be  easily  solved,  and  their 


solutions  are 


u(t)  = XQCOsh(kt)  - cyQSinh(kt) 
v(t)  = y^cosh(kt)  - c~^x  sinh(kt) 


(2.11) 


with  c = (a/b)^^^  and  k = (ab)^^^. 


'aii(t)  \ /a  a ® \ [ d^exp(kt)  + d2exp(-kt)  + d^exp(-2kt) 

0i2(t)  = — k -k  0 dj^exp(2kt)  + d^exp(kt)  + dgexp(-kt) 

8k 

yb  b -b  y ycL^exp(kt)  + d0exp(-kt)  + d^ 

with 

di  = j(yo(l-c)  + Xq(1-c'^)),  dg  = yo(l+c)  + Xq(1+c"^), 

= -^(yo(^-2c)  + Xq(4-2c"^)),  d^  = |(yQ(4-2c)  + Xq(4-2c‘^)), 

d^  = -(yQ(l-c)  + Xq(1-c“^)),  dg  = -|(yQ(l+c)  + Xq(1+c“^)), 

(Lj.  = 2(yQ(l+c)  - Xq(1+c'^)),  dg  = -2(yQ(l-c)  - Xq(1-c‘^)), 


<^9  = -^(cyQ  - c-1xq) 
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The  diffusion  approximation  thus  gives  the  time  evolution  of 
(X(t),  Y(t))  as  the  sum  of  the  deterministic  process  given  by 
(2.11)  and  a stochastic  noise  process  which  has  been  completely 
characterized.  This  illustrates  that  the  deterministic  approach 
represented  by  (2.1)  is  a first-order  approximation  to  the  Lan- 
chester  process.  A second-order  approximation  is  obtained  by 
adding  a stochastic  noise  process.  The  diffusion  approximation 
methodology  provides  a link  between  the  deterministic  and  sto- 
chastic models. 

The  reader  ^ould  be  aware  that  the  importance  of  using  a 
stochastic  analysis  will  depend  on  the  size  of  N.  For  very  large 
values  of  N,  the  deterministic  term  will  dominate  the  stochastic 
term.  In  such  cases,  the  deterministic  approach  will  be  completely 
adequate.  On  the  other  hand,  when  N is  small,  say  less  than 
500,  wide  fluxuations  are  possible.  To  illustrate,  in  a square 
law  case  with  a=b=  .05,  ~ ^0  ^ 300,  after  15  time  units 

the  X force  level  could  be  as  low  as  102  or  as  high  as  I85. 

If  Xq  = yQ  = 50,  the  X force  level  may  be  between  1 and  27* 
The  deterministic  model  gives,  therefore,  very  misleading  results. 

The  reader  should  be  aware  of  an  important  limitation  of  the 
diffusion  approximation  approach,  the  problem  of  boundaries.  One 
must  be  aware  that  the  (X(t),  Y(t))  process  is  confined  to  the 
first  quadrant  with  the  x and  y axes  being  absorbing  barriers. 
The  bivariate  normal  distribution  will  be  positive  over  the  entire 
plane,  while  the  process  can  never  leave  the  first  quadrant.  The 
approximation  will  be  trustworthy  only  if  most  of  the  mass  of 
the  normal  distribution  is  in  the  first  quadrant. 
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3.  An  Analysis  of  Specific  Models 

In  this  section,  we  illustrate  the  method  of  diffusion 
approximations  by  applying  it  to  a number  of  different  Lanchester 
processes.  The  models  chosen  for  study  have  been  selected  because 
they  have  classical  significance  (the  linear  law),  illustrate 
mathematical  points  (Helmbold's  general  attrition  structure),  or 
involve  complex  modelling  (the  heterogeneous  and  engagement  models 
of  Karr).  The  reader  will  see  how  diffusion  approximations  pro- 
vide a unifying  tool  for  the  study  of  stochastic  Lanchester  pro- 
cesses. 

A.  The  Linear  Law 

The  linear  law  arises  when  the  functions  f and  g in 
(2.2)  are  f(i,j)  = aiJ/N  and  g(i,j)  = biJ/N.  Equation  (2.6) 
becomes 

dUjj(t)  = -N^/^(u'(t)  + au(t)v(t))dt  + 0(N"^/^)dt 

-(av(t)Ujj(t)  + au(t)Vjj(t)  )dt 

+ (au(t)v(t)  + 0(n"^^^) )^/^dW.  (t) 

^ (3.1) 

= -N^^^(v'(t)  + bu(t)v(t))dt  + 0(N"^/^)dt 
-(bv(t)Ujj(t)  + bu(t)Vjj(t)  )dt 
+ (bu(t)v(t)  + 0(N‘^^^) )^/^dW2(t) 


One  may  again  appeal  to  Theorem  k of  Barbour  (1974)  to 
prove  that  if 


u' ( t)  « -au( t) v( t) 
v'  ( t)  * -bu( t) v( t) 


(3.2) 

u(0)  = X(0)/N,  v(0)  = Y(0)/N, 
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then  [(Ujj(t),  Vjj(t)),  t>0}  converges  weakly  to  {(U(t),  V(t)), 
t>0]  which  satisfies  (2.7)  with 


At  = 


-av(t)  -au(t) 
-bv(t)  -bu(t)  I 


and 


5t  = 


(au(t)v(t) ) 


1/2 


(bu(t)v(t)) 


1/2  • 


The  diffusion  approximation  gives  (X(t),  Y(t))  w N(u(t),  v(t)), 

+ N^^^(U(t),  V(t))  where  u(t)  and  v(t)  satisfy  the  usual  Lanches- 
ter  equations  (5.2).  The  distribution  of  (X(t),V(t))  will  thus  be 
approximately  bivariate  normal  with  the  usual  deterministic  result  as  a 
mean  and  covariance  where  can  be  derived  from  (2.9)  or 

equivalently  (2.10)  with 


1 -2av( t) 

-2au( t) 

0 

f au(t)v(t) 

^t  = 

' -bv(t) 

-(av(t)  + bu(t) ) 

-au(t)  , = 

1 

0 

1 ° 

-2bv( t) 

-2bu( t) 

1 bu( t) v( t) 

The  deterministic  equations  (5.2)  can  be  solved  in  closed 
form;  however,  the  covariance  equations,  having  time  dependent 
coefficients,  will  typically  not  be.  Nevertheless,  (5*5)  can  be 
solved  by  straightfoT^/^ard  numerical  techniques.  If  numerical 
techniques  are  required,  then  the  f and  g functions  can  be 
made  to  depend  explicitly  on  time.  Alternative  expressions  for 
the  solution  of  (5.5)  are  given  in  [1],  Chapter  8. 

It  is  important  to  note  that  using  diffusion  approximation 
methodology  reduces  the  problem  of  studying  stochastic  Lanchester 
processes  to  that  of  solving  ordinary  differential  equations  which 
arise  from  the  deterministic  approach  and  the  covariance  matrix. 

As  such,  the  stochastic  approach  is  In  principle  no  harder  than 
the  deterministic  approach. 
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B.  The  General  Attrition  Strxicture  of  Helmbold 

In  a paper  by  Helmbold  [8],  a general  attrition  structure 
was  proposed  which  incorporates  inefficiencies  of  scale  associated 
with  a vary  large  force  fighting  a relatively  smaller  one.  The 
attrition  structure  proposed  consists  of  replacing  (2.1)  by 

x'(t)  = -ax^(t)y^~'^(t),  x(0)  = x^,  y(0)  = y^ 

y'(t)  = -by'^(t)x^''^(t),  (5.^) 

and  0 < c < 1 . 


The  associated  stochastic  version  consists  of  transitions 
and  transition  rates  similar  to  (2.2)  given  by 


t t + dt 

(i,J)  ^ (i-l,J) 

(i,J-l) 
(i,J) 

all  others 


Rate 

ai'^J^-'^dt  + o(dt) 
bj^i^"'^dt  + o(dt) 

1 - (ai'^J^'*^  + bi^i^"'^)dt  + o(dt) 


o(dt) . 


The  non-linear  form  of  the  transition  rates  precludes  any 
possibility  of  an  exact  expression  for  the  distribution  of  X(t), 
Y(t))  being  derived.  Nevertheless,  the  process  can  be  approximated 
by  using  s.d.e.'s  which  are  obtained  by  matching  the  first  and 
second  moments.  Specifically,  we  write 

dX(t)  = -a(X(t))‘^(Y(t))^“^dt  + (a(X(t))'^(Y(t))^'^)^/^dWi(t) 
dY(t)  - -b(X(t))^'^(X(t))^dt  + (b(X(t))^‘^(Y(t))^)^/^dW2(t), 

with  X(0)  = Xq,  Y(0)  = Yq  a.s. 

We  once  again  introduce  Uj^(t)  and  Vj^(l)  defined  in 
section  2 and  using  Ito's  Lemma  write 
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dUj^(t)  = -N^/^(u'(t)  + aN‘^/^(Nu(t)  +N^/^Uj^(t))'^(Nv(t)  +N^/^Vj^(t))^“^)dt 
+ (aN"^(Nu(t)  +N^/^Uj^(t))^(Nv(t)  + t)  )^"^)^/^dW^(  t) 

dVj^(T)  = -N^/^{V  (t)  +bN'^/^(Nu(t)  +N^/^Uj^(t))^"^(Nv(t)  + 

+ (bN"^(Nu(t)  +N^/^Ujj(t))^“'^(Nv(t)  +N^/^Vjj(t))^)^/2dW2(t) 

(3-7) 

The  factors  such  as  (Nu(  t)  + t) ) ^ can  be  expanded 

in  a power  series.  Upon  performing  these  expansions,  collecting 

-1/2 

terras,  and  dropping  all  terms  which  are  0(N~  ' ) we  find 

/dUj^(t)\  _ ^/2 


dV.(t) 


v'(t)  + bu  (t)v  (t) 


ca(Ziii) 

u(t) 

(l.c)b(lIlX)' 

u(t) 


a{l-o)(ii5i) 

u(t) 


,1-c- ,.1/2 


(au‘"(t)v^'^(t)) 


cb(:^) 

u(t) 


(bu^"^(t)v^(t))- 


Uij(t) 

v^(t) 


(3.8) 


dW^(t) 

dW2(t) 


If  u and  V satisfy 

u'(t)  - -au‘^(t)v^"^(t)  u(0)  = Xq/N 

v'(t)  = -bu^"'^(t)v^(t)  v(0)  = y^/N 

then  t(Ujj(t),  Vjj(t)),  t>0]  converges  weakly  to  ((U(t),  V(t)), 
t > 0}  which  satisfies  (2.7)  with 
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= 


St  = 


u(t) 

(l-c)b(^)' 
V u(t) 


a(l-c)(^)’' 

u(t) 

cb(lLD.)''“^ 

u(t) 


(au^(t)v^"'^(t))^/2 


(bu^"‘=(t)v'^(t))^/2 


The  mean  of  the  stochastic  process  is  given  by  N(u(t), 
v(t)),  while  the  covariance  is  given  by  where  satisfies 

(2.9)  or  equivalently  (2.10)  with 


2car^~*^  ( t) 


2a(l-c) r“^(  t) 


= -l(l-c)br^(t)  c(ar^“'^(t)  + br^“^(t) ) 

\ 0 2{l-c)br'^(t) 


a(l-c)r'^(t) 
2cbr^"^(t)  i 


with  r(t)  = v(t)/u(t)  and 


= (au^(t)v^"^(t),  0,  bu^"^(t)v*^(t)  )^. 


Again,  (2.10)  must  be  solved  numerically,  but  the  stochastic 
process  problem  has  been  reduced  to  solving  two  systems  of  ordi- 
nary differential  equations,  even  in  the  case  of  highly  nonlinear 
transition  rates. 


C.  Homogeneous  Linear  Law  Process  With  Engagements 

A number  of  interesting  complex  stochastic  models  have  been 
introduced  by  Karr  [10].  The  first  which  we  consider  is  that  of 
a linear  law  process  with  engagements.  In  this  process,  each 
side  has  a homogeneous  force.  The  two  sides  suffer  attrition  by 
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one-on-one  engagements  or  encounters.  This  means  that  each  side 
can  be  thought  of  as  being  divided  into  two  parts  - a part  which 
is  involved  in  engagements  and  a part  seeking  such  involvements. 
To  make  the  stochastic  process  Markovian,  one  must  assume  that 
engagements  last  for  an  exponential  length  of  time  with  parameter 
U.  Furthermore,  each  engagement  has  three  possible  outcomes: 
side  i wins  with  probability  p^^,  i=l,2  or  both  survive  with 
probability  p^  (p^ + p^ + p^ = 1) . Engagements  are  assumed  to 
develop  at  a rate  proportional  to  the  product  of  the  size  of  the 
two  unengaged  force  levels.  In  this  regard,  the  model  has  a 
certain  similarity  to  the  stochastic  epidemic  model. 

There  are  three  classifications;  X(t)  and  Y(t),  the 
unoccupied  force  levels  at  time  t,  and  Z(t),  the  engagements 
at  time  t.  The  [(X(t),  Y(t),  Z(t)),  t>0]  process  can  be 
modelled  as  a trivariate  Markov  process  with  transitions  and 
rates  given  by 


t t -I-  dt 

(i,J,k)  ^ (i+l,J,k-l) 

(i,J+l,k-l) 
(i+l,J+i,k-l) 
(i-l,J-l,k+l) 
all  others 


Rate 

P2(ikdt  + o(dt) 

p^tikdt  + o(dt) 

p^ukdt  + o(dt)  (3.10) 

\N"^ijdt  + o(dt) 
o(  dt) 


where  N = X(0)  + Y(0)  + Z(0)  and  is  assumed  to  be  large.  The 
approximate  s.d.e.'s  can  be  written 
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dX(t)  = (P2  + P3VZ(t)dt  - \N“^X(t)Y(t)dt  + (p^uZC  t ) t ) 

+ (P3UZ(t))^/^dW2(t)  - (XN"^(t)Y(t))^/^dW^(t) 

dY(t)  = (p^ + Pj)nZ(t)dt  - XN‘^X(t)Y(t)dt  + (p^uZ(t))^/^dW2^(t) 

+ (P3UZ(t))^^^dW2(t)  - (\N"^X(t)Y(t))^/^dW^(t) 

dZ(t)  = -uZ(t)dt  + \N“^A(t)Y(t)dt  - (p^uZ(t))^/^dW|^(t) 

- (P2uZ(t))^/^dW^(t)  - (p^uZ(t))^/2dW^(t) 

+ (xN"^X(t)Y(t))^/^dW^(t).  (3.11) 


Once  again,  we  introduce  the  normalized  process  {(Uj^(t), 
Vj^(t),  Wjj(t)),  t>0}  where  Ujj(t)  = t ) - Nu(  t ) ) , V^{t)  = 

jj-l/2(Y(t)  . Nv(t)),  and  Wj^(t)  = Zj^(  t ) - Nw(  t ) ) . This 

process  will  be  a trivariate  diffusion  with  s.d.e.  which  can  be 
derived  using  Ito' s Lemma.  If  one  lets  N oo  and  assumes  u,  v, 
and  w satisfy 


u'(t)  = (P2  + P3)uw(t)  - \u(t)v(t) 

v'(t)  = (p^ +P^)|iw(t)  - xu(t)v(t)  (3-12) 

w'(t)  = -uw(t)  + Xu(t)v(t), 


then  [(Ujj(t),  Vjj(t),  Wjj(t)),  t>0}  converges  weakly  to  [(U(t), 
V(t),  W(t)),  t>0}  which  satisfies 


dU(t) 

dV(t) 

dW(t) 

I 


-Xv(t)  -Xu(t)  (P2  + P3)uw(t) 
-Xv(t)  -xu(t)  (Pi+P^)uw(t) 
Xv(t)  Xu(t)  -uw(t) 


U(t)^ 

V(t) 

W(t) 


dt 


( continued) 


(3.13) 


20 


0 

(P^hw(t) 

-U 

u(t)v(t))^/2  ^ 

0 

(p^^w(t))^/2 

(p^^iw(t)  )^/^ 

-(X 

u( t) v( t) 

-(P24w(t))^/^ 

-(PlUw(t))^/^ 

-(Pj^w(t))^/^ 

(X 

u(t)v(t))^/^  ^ 

or  in  matrix  form 

dU(t)  = A^(t)dt  + B^dW(t)  with  U(0)  = 0. 

The  limiting  stochastic  process  is  linear,  so  the  process  will 
have  trivariate  normal  marginals.  The  mean  of  U(t)  is  0. 

The  covariance  at  time  t,  is  the  unique  non-negative 

definite  solution  of  (2.9) • Both  equations  (J.12),  the  classi- 
cal Lanchester  equations  for  this  model,  and  (2.9)  must  be 
solved  numerically.  Since  a numerical  solution  is  required,  the 
problem  is  essentially  no  harder  to  solve  if  X,  u,  and  the 
p' s are  made  explicitly  time  dependent. 

The  diffusion  approximation  thus  characterizes  the  attri- 
tion process  as  being  Gaussian  with  specified  mean  and  covari- 
ance fiinctlon,  namely  N(u(t),  v(t),  w(t))  and  NZ^.  The 
multivariate  case  thus  poses  no  additional  problems  for  this 
approximate  methodology  whereas  the  forward  equations  become 
more  difficult  to  solve. 


D. 


Heterogeneous  Linear  Law 

Karr  [10]  has  introduced  a variety  of  different  heterogeneous 
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models.  These  are  characterized  by  each  side  having  a number 
of  different  types  of  forces.  Each  force  type  on  each  side 
engages  each  force  type  of  the  other  side.  These  engagements 
can  be  of  any  type;  linear,  quadratic,  etc.  Thus  many  different 
attrition  processes  are  going  on  simultaneously.  A rich  variety 
of  different  models  are  possible,  and  we  illustrate  the  diffusion 
approximation  analysis  in  a few  special  cases.  We  present  only 
the  case  where  both  sides  have  two  distinct  types  of  forces. 

The  analysis  can  be  easily  generalized  to  the  case  where  each 
side  has  many  types.  The  reader  should  be  aware,  however,  that 
while  the  mathematics  is  straightforward,  the  diffusion  approxi- 
mation approach  is  accurate  only  when  the  number  of  each  type 
is  reasonably  large.  As  the  number  of  types  increases,  this 
assumption  will  become  doubtful. 

We  consider  two  force  level  process  X^(t)  and  X2(t)  for 
one  side  and  Y^(t),  Y2(t)  for  the  other  at  time  t.  For  a 
completely  linear  case,  one  would  treat  [(X^(t),  X2(t),  Y^(t), 
Y2(t)),  t>0}  as  a continuous  time  Markov  chain  with  transi- 
tions and  transition  rates  given  by 

t t + dt 

(i,J,k,l)  ^ (i-l,J,k,l) 

(i,J-l,k,l) 

(i,J,k-l,l) 

(i,J,k,l-l) 

The  diffusion  approximation  analysis  would  begin  by 
approximating  the  process  using  s.d.e.'s  given  by 


Rate 

i(c^k+  C2l)dt  + o(dt) 

J(d.,k+ dol)dt  + o(dt) 

12'  (J.lit) 

k(e^i + 62J)dt  + o(dt) 
l(f^i + f2J)dt  + o(dt). 
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dX^(t)  = -X^(t)(c^N"^y^(t)  4-  C2N“^Y2(t))dt 
dX2(t)  = -X2(t)(d^N“^Y^(t)  + d2N"^Y2(t))dt 

+ (d^N'42(t)Y^(t))^/^dW^(t)  + (d2N"^X2(t)Y2(t))^/2^y^(^) 
dY^(t)  = -Y^(t)(e^N"^X^(t)  + e2N"^X2(t))dt 

+ (e^N'\(t)Y^(t))^/^dW^(t)  + (e2N"^X2(t)Y^(t))^/2^^j^(^^ 

dY2(t)  = -Y2(t)(f^N'^X^(t)  + f2N"^X2(t))dt 

+ (f^r^X^(t)Y2(t))^/^dW^(t)  + (f2N"42(t)Y2(t))^/2^^^(^j 

where  N = X^(0)  + X2(0)  + Y^(0)  + Y2(0).  (5.15) 

The  analysis  proceeds  in  the  same  fashion  as  before.  The 
normalized  diffusion  processes  “ (X^(t)  - Nu^( t) 

and  = (Yj_(t)  - NVj^(  t)  i=l,2  are  introduced  and 

the  s.d.e.  governing  their  behavior  is  derived  using  Ito's 
Lemma.  Finally,  one  lets  N-^oo  assuming  u^  and  v^,  i=l,2 
satisfy  the  ordinary  differential  equations 

uj_(t)  = -u^(t)(c^v^(t)  + c^v^it)) 
u^(t)  = -U2(t)(d^v^(t)  + d2V2(t)) 

(5.16) 

vj_(t)  = -v^(t)(e^u^(t)  + e2U2(t)) 
v^(t)  = -V2(t)(f^u^(t)  + f2U2Ct)). 

The  t(U^jj(t),  U2jj(t),  V^N(t).  V2j,(t)),  t>0}  process 
converges  weakly  to  {(U^(t),  U2(t),  V^(t),  Vg(t)),  t>0)  which 


satisfies 
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where  the  t dependence  of  and  has  been  supressed. 

The  analysis  of  (^.17)  is  straightforward,  since  it  is  linear. 
The  distribution  at  time  t will  be  multivariate  normal  with  mean 
0 and  covariance  specified  by  (2.9).  It  follows  that  (X^(t), 

X2(t),  Y^(t),  Y2(t))  will  have  approximately  a normal  distribu- 
tion with  mean  N(u^(t),  U2(t),  v^(t),  V2(t))  and  covariance 
N2^.  It  appears  that  (^5.16)  and  (2.9)  must  be  solved  numerically. 

E.  Heterogeneous  Quadratic  Law 

The  analysis  presented  in  can  be  easily  modified  to  give 
an  analysis  of  the  heterogeneous  quadratic  case.  Transitions 
given  by  are  replaced  by 


t t + dt 

(i,J,k,l)  ^ (i-l,J,k,l) 

(i,  j-l,k,l) 
(i,j,k-l,l) 
(i, j,k,l-l) 


rate 

(c^k  + c^ljdt  + o(dt) 

(d.k  + d..l)dt  +o(dt) 

^ ^ (3.18) 

(e^i  + 62J)dt  + o (dt) 

(f^i  + f2J  )<lt  + o (dt) 


The  deterministic  equations  given  by  (3.16)  are  replaced  by 


u^(t)  = -(c^v^(t)  + C2V2(t)) 
u^(t)  = -(d^v^(t)  + d2V2(t) ) 
vj_(t)  = -(e^u^(t)  + e2Ug(t) ) 
v^(t)  = -(f^u^(t)  + f2U2(t)). 


(3.19) 


These  equations  can  be  easily  solved  using  standard  methods. 
The  functions  u^,  Ug,  v^,  and  can  be  written  as  a sum  of 

exponentials. 

The  stochastic  noise  term  will  satisfy  a s.d.e.  similar  to 
(3.17).  It  will  be  given  by 
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where  the  t dependence  of  and  v^,  i=l,2  has  been 

supressed. 

The  mean  of  (Uj^(t),  U^it) , V^(t),  V^it))  is  0 while 
the  covariance  matrix  is  the  solution  of  (2.9).  This 

matrix  can  be  calculated  in  closed  form  and  will  be  a mixture 
of  exponentials. 

The  analysis  presented  in  3D  and  3E  can  be  easily  generalized. 
Mixtures  of  attrition  structures  or  indeed  far  more  complicated 
structures  can  be  easily  accomodated.  The  latter  two  sections 
are  meant  only  to  illustrate  the  generalization  to  heterogeneous 
structures. 

4.  Numerical  Results 

Extensive  numerical  studies  were  carried  out  to  assess  the 
accuracy  of  the  diffusion  approximation  methodology.  The  method 
is  based  on  an  asymptotic  analysis  where  X(0)  and  Y(0)  are 
assumed  to  be  proportional  to  N and  N^«.  Even  so,  we  have 
claimed  that  the  method  provides  good  results  for  values  of 
X(0)  and  Y(0)  as  small  as  30.  We  show  only  two  situations 
in  this  section,  both  based  on  the  "quadratic  law";  however,  they 
are  representative  of  the  many  cases  considered. 

A set  of  parameter  values  and  starting  values  was  chosen  as 
well  as  a time  T.  The  distribution  is  predicted  by  the  diffusion 
approximation  methodology  to  be  bivariate  normal  with  specified 
mean  and  covariance  structure.  To  assess  this  prediction,  the 
actual  Markov  chain  given  by  (2.2)  was  simulated  up  to  time  T 
6000  times.  The  simulation  results  were  compared  with  the 
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approximation  in  several  ways.  First,  the  means,  variances  and 
covariances  were  compared.  Second,  the  simulated  distributions 
of  X(T)  and  Y(T)  were  plotted  on  normal  probability  paper 
to  determine  the  accuracy  of  the  normal  distribution  approxima- 
tion. Third,  the  (X(T),  Y(T))  pairs  were  linearly  transformed 
and  the  components  plotted  on  normal  probability  paper  to  deter- 
mine whether  the  joint  distribution  was  really  bivariate  normal. 


a = .05, 

CT 

II 

o 

II 

15 

(xo,yo) 

X(T)  - mean 

Y(T) 

- mean 

2 

D 

|S-D|/S 

S 

D 

IS-Di/S 

(20,20) 

9.48 

( .041) 

9.^5 

.0032 

9.4o 

( .002) 

9.45 

.0053 

(25,25) 

11.85 

(.069) 

11.81 

.0034 

11.87 

(.061) 

11.81 

.0051 

(30,30) 

14.21 

(.041) 

14.17 

.0028 

14.16 

(.011) 

14.17 

.0007 

(40,40) 

18.98 

(.190) 

18.89 

.0047 

18.93 

(.069) 

18.89 

.0021 

(50,50) 

23.64 

(.118) 

23.62 

.0008 

23.58 

(.251) 

23.62 

.0017 

X(T)  - Standard  Deviation  Y(T)  - Standard  Deviation 


X 

0 

0 

S 

D 

1 S-Dl/S 

s 

D 

IS-Dl/S 

(20,20) 

3.87 

(.053) 

3.98 

.0284 

3.9^ 

( .026) 

3.98 

.0102 

(25,25) 

4.42 

(.050) 

4.45 

.0068 

4.39 

(.064) 

4.  45 

.0137 

(30,30) 

4.84 

(.004) 

4.88 

.0083 

4.83 

(.091) 

4.88 

.0104 

(40,40) 

5. 66 

(.136) 

5.63 

.0053 

5.  64 

(.027) 

5.63 

.0018 

(50,50) 

6.27 

(.031) 

6.30 

.0048 

6.31 

(.033) 

6.30 

.0016 

(X(T),  Y(T))  Correlation 

(xo,yo) 

S 

D J 

S-Dl/S 

(20,20)  -.69^ 

(.005) 

-.686 

.0115 

(25,25)  -.684 

(.015) 

-.686 

.0029 

(30,30)  -.683 

( .001) 

-.686 

.0044 

(4o,4o)  -.689 

( .006) 

VO 

00 

VO 

.0044 

(50,50)  -.687 

(.009) 

-.686 

.0015 
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a = .075,  b = .050,  T = 10 
X(T)  - mean 

(Xo,yo)  S D 

Y(T)  • 

|S-Di/S  S 

- mean 
D 

1 S-Dl/S 

(50,20) 

40.18 (.003) 

4o.l6 

.0005  6.78(.056) 

6.72 

. 0088 

(75,30) 

6o.32( .029) 

60.24 

.0013  10. 08 (.026) 

10.09 

.0010 

(100,40) 

8o.34(.o66) 

80.32 

.0002  13.45(.155) 

13.45 

.0000 

(125,50) 

100.32(.030) 

100.41 

.0009  16.86(.o89) 

16.81 

.0030 

(250,100) 

200.96( .149) 

200.81 

.0007  33.70(.008) 

33.62 

.0024 

X(T)  - Standard  Deviation 

Y(T)  - standard  Deviation 

X 

o 

o 

S E 

!_  is 

-D|/S  S D 

1 S-D 

l/s 

(50,20) 

3.59(.021)  3. 

68 

0251  3.58  (.035)  3.84 

.0726 

(75,30) 

4.43(.055)  4. 

51  . 

0181  4.61  (.074)  4.71 

.0217 

(100,40) 

5.26(.o49)  5. 

21 

0095  5.40  (.096)  5.43 

.0056 

(125,50) 

5.8o(.070)  5. 

82  . 

0034  6.08  (.088)  6.08 

.0000 

(250,100) 

8.25(.238)  8. 

24 

0012  8.62  (.054)  8.59 

.0035 

(X(T),  Y(T))  Correlation 

( 

S 

D 

IS-Dl/S 

(50,20) 

-.512(.011)  - 

.526 

.0273 

(75,50) 

-.534(.008)  - 

.526 

.0150 

(100,40) 

-•547(.023)  - 

.526 

.0384 

(125,50) 

-•518(.038)  - 

.526 

.0154 

(250,100) 

-.527(.006)  - 

.526 

.0019 

TABLE  2 


29 

The  columns  marked  S and  D in  Tables  1 and  2 represent 
the  simulation  and  diffusion  approximation  results  respectively. 
The  quantity  |S-D|/S  gives  the  percentage  error.  In  the  S 
column,  for  each  starting  value  an  estimate  and  its  standard 
error  are  presented.  Often  the  quantity  [S-D|  is  smaller  than 
the  standard  error. 

Tables  1 and  2 illustrate  the  extraordinary  accuracy  of 
the  diffusion  approximation  method  even  for  starting  values  as 
small  as  20.  The  mesuis  are  off  by  no  more  than  1/2  percent 
error,  a quantity  easily  explained  by  the  standard  error  of 
the  simulation.  Similarly  an  error  of  less  than  5^  is  noted 
for  the  standard  deviations  and  correlations  and  this  drops  to 
below  1/2  percent  as  the  starting  values  increase.  The  diffu- 
sion approximation  methodology  predicts  moments  very  accurately. 

Figures  1,  2,  and  ^ give  normal  probability  plots  of  X(T) 
and  Y(T)  for  various  parameter  settings.  The  plots  are  shown 
for  cases  with  a small  starting  values  ((20,20)  or  (50,20)) 
and  with  a moderate  starting  values  ((50,50)  or  (125,50)). 

The  plots  for  small  starting  values  exhibit  some  curvature 
indicative  of  non-normalityj  however,  it  is  rather  slight.  This 
means  that  probabilities  calculated  from  a normal  distribution 
will  be  reasonably  accurate. 

The  cases  of  moderate  starting  values  show  very  strong 
normality  (making  allowances  for  sampling  variability).  These 
plots  clearly  indicate  that  when  Xq  and  y^  are  50  or  more 
excellent  accuracy  will  be  achieved  using  the  diffusion  approxi- 


mation. 
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The  plots  of  the  transformed  X(T)  and  Y(T)  variables 
are  not  shown  but  gave  results  in  accord  with  those  illustrated 
in  Figures  1,  2,  and  3* 

The  figures  also  sejrve  to  illustrate  the  importance  of 
using  a stochastic  analysis  rather  than  a deterministic  analysis. 
In  Figure  1 with  (x^,  Yq)  = (20,20),  a deterministic  analysis 
gives  (x(l5),  y(15))  = (9-^5,  9-^5)-  The  graph  shows  P(X(15) 
<6)  = .22  and  P(  X(15)  > 12 ) = . 22,  so  P(T  < X(15)  > H ) = . 56. 

The  deterministic  value  9*^5  ls»  therefore,  quite  misleading 
in  view  of  the  large  variability.  The  case  with  (Xq,  y^)  = 
(50,50)  is  even  more  pronounced.  Here  P(21 < X(T)  < 29)  = . 52, 
while  the  deterministic  value  is  23.62. 

It  is,  therefore,  clear  that  unless  (xq,  y^)  is  very 
large  and  the  attrition  essentially  deterministic,  a stochastic 
analysis  must  be  employed.  In  view  of  the  accuracy  and  con- 
venience of  the  diffusion  approximation  it  appears  to  be  an 


ideal  method. 


Figure 
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